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Studying time-dependent behavior in lasers is analytically difficult due to the saturating non¬ 
linearity inherent in the Maxwell-Bloch equations and numerically demanding because of the 
computational resources needed to discretize both time and space in conventional FDTD approaches. 
We describe here an efficient spectral method to overcome these shortcomings in complex lasers 
of arbitrary shape, gain medium distribution, and pumping profile. We apply this approach to 
a quasi-degenerate two-mode laser in different dynamical regimes and compare the results in the 
long-time limit to the Steady State Ab Initio Laser Theory (SALT), which is also built on a spectral 
method but makes a more specific ansatz about the long-time dynamical evolution of the semiclassical 
laser equations. Analyzing a parameter regime outside the known domain of validity of the stationary 
inversion approximation, we find that for only a narrow regime of pump powers the inversion 
is not stationary, and that this, as pump power is further increased, triggers a synchronization 
transition upon which the inversion becomes stationary again. We provide a detailed analysis of 
mode synchronization (aka cooperative frequency locking), revealing interesting dynamical features 
of such a laser system in the vicinity of the synchronization threshold. 


Lasers are very rich dynamical systems which exhibit 
various time-dependent phenomena characteristic of non¬ 
linear systems such as phase and mode locking, self¬ 
pulsing and breathing, and generally, spatio-temporal 
pattern formation and dynamical chaos. Almost all these 
effects can be understood and quantitatively studied using 
the semiclassical laser theory in the form of Maxwell-Bloch 
(MB) equations [IHS], a set of coupled non-linear equations 
for the space- and time-dependent electric field amplitude 
E(r,t), and the polarization and inversion of the gain 
medium P{r,t) and D(r,t). Early work made abundant 
use of spectral methods, where the field amplitudes en¬ 
tering the MB equations are expanded in a complete 
basis of spatial modes, reducing MB equations to a set 
of coupled non-linear ordinary differential equations for 
time-dependent amplitudes. These early theoretical inves¬ 
tigations made a number of simplifying assumptions on 
the spatial aspects of the problem. The lasing modes were 
assumed to be simple (uniform, trigonometric, or gaus- 
sian) and unmodified from their passive cavity modes, 
and the openness (optical leakage) was taken into ac¬ 
count phenomenologically. While these assumptions are 
sufficiently general to reproduce qualitatively almost all 
features of laser dynamics in macroscopic cavities, new 
laser systems have emerged in the past two decades that 
raised questions not easily addressable by these spectral 
approaches. 

Most novel laser systems are motivated by their deploy¬ 
ment as compact and tunable light-sources for on-chip 
applications [3]. Typically, these lasers feature complex 
sub-wavelength patterning of the cavity volume to employ 
light-confinement mechanisms that are based on optical 
interference (photonic band gap materials that may or 
may not include optical defects, random lasers) and/or 
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total internal reflection (whispering gallery lasers, wave- 
chaotic lasers). Therefore these lasers feature spatially 
complex modes, typically in highly open geometries. In 
some cases, such as weakly scattering random lasers, it is 
not even clear where the boundary of the system is, and 
even what a mode means [5] . In addition, many solid-state 
lasers are subject to spatially non-uniform pumping con¬ 
ditions and feature strong modal interactions [SHin]. All 
these conditions can be modeled by appropriately setting 
up the original Maxwell-Bloch equations and solving the 
resulting non-linear partial differential equations (PDEs) 
in time-domain through various finite-difference-based nu¬ 
merical methods HIHlsl. A number of such powerful com¬ 
putational methods have been developed and employed 
to investigate the dynamics of complex laser systems, ei¬ 
ther solving the full set of MB equations [I4lfl6] , or the 
parabolic version obtained upon a slowly varying envelope 
(SVE) approximation in the time-domain, the so-called 
Schrodinger-Bloch (SB) equations [17] . 

A more recent approach, the Steady-state Ab-initio 
Laser Theory (SALT) [181 [19], overcomes the rather ex¬ 
pensive discretization of the spatial domain of a complex 
laser system in MB/SB-FDTD solvers by taking a spectral 
approach. The field amplitudes E{r,t) are expressed in 
the Constant-Flux (CF) basis [18], a set of non-Hermitian 
modes that exactly describe the steady-state field distribu¬ 
tion in a finite and open domain under harmonic driving 
conditions |2()| . There are a number of advantages pro¬ 
vided by this approach. (1) The steady-state multi-mode 
solution (to be defined precisely below) in the asymptotic 
infinite-time limit is obtained directly, without resorting 
to a time-domain simulation, (2) The exact solution of the 
MB equations in the steady-state is obtained through a 
modular two-stage procedure: in the first stage the linear 
problem corresponding to the determination of a CF basis 
is solved, and in the second stage this information is used 
to solve a set of algebraic transcendental equations [20]. 
This allows the separation of spatial complexity (handled 
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as a linear problem) from the computational non-linear 
problem, and perhaps more importantly obviates the need 
for the computational implementation of boundary condi¬ 
tions through various PML-variety approaches [5T]. (3) 
SALT is also flexible enough to effectively account for spa¬ 
tially non-uniform pumping conditions [IHl [El 1221 ESj , (4) 
can directly provide the farfield electric field distribution 
and spectrum m, and (5) over the past decade provided 
unique semi-analytic insight to fundamental problems 
in laser physics [El lEl EM] that is harder to attain 
through brute-force computational approaches. 

Despite the success of SALT in the treatment of com¬ 
plex laser systems, it has certain well-known limitations. 
The key assumption of the theory is the stationarity of 
the inversion D{r,t) [17] (or more generally, the level 
populations |28jL The inversion is however never truly 
stationary, but in a certain regime of parameters the non¬ 
stationary corrections are systematically very small and 
can be neglected. To be more specific, the non-stationary 
corrections as discussed in detail below, are order 
7 ||/A where A is the smallest frequency difference of the 
lasing modes (typically slightly different from the free 
spectral range of the cold cavity) and yy is the inver¬ 
sion relaxation rate. However, A depends on the pump 
strength, and at larger powers can become smaller than 
7 || due to non-linear effects. As a consequence, the cor¬ 
rections to SALT are not going to be small, and can lead 
to qualitatively different behavior. As shown here, such 
a scenario can take place under unusual circumstances 
where the spectrum of the cavity contains quasi-doublets 
(typically protected through a discrete spatial symmetry 
of the cavity) that are spectrally spaced apart at a dis¬ 
tance (A 2 ) that is larger than the splitting of the doublets 
(A), as shown in Fig. Under such circumstances the 
lasing modes of the doublet-pair most favored by the gain 
curve can lock to each other and synchronize as pump 
power is increased through an effect called ’cooperative 
frequency locking’ |53|. Yet even then, as we will show, 
the Stationary Inversion Approximation (SIA) fails only 
in a very limited pump power range near the synchroniza¬ 
tion threshold, and is valid for most of the pump power 
range below and above this threshold. 

Thus it is of interest not only to understand the valid¬ 
ity of the SIA under various circumstances, but also to 
develop a spectral method that is in principle not limited 
by any approximations such as the SIA. Such a technique 
should be able to capture any intrinsically dynamical be¬ 
havior of complex lasers. We present such a technique 
in this article, and discuss precisely how SIA and thus 
SALT may fail in certain limited parameter regimes. 

Just as SALT, the new Constant Flux Time Domain 
(CFTD) technique presented here provides a versatile 
tool for calculating lasing thresholds, spectra, and modal 
distributions in the multi-mode regime for complex lasers 
including random (El) semiconductor [30], photonic crys¬ 
tal surface-emitting [31] , and photonic molecule lasers [32] . 
Unlike SALT, it can capture transient regimes, locking 
and synchronization, various dynamical instabilities [33] . 


as well as dynamical chaos and generally, spatio-temporal 
pattern formation. 

In Section I we provide an overview of our theoretical 
approach, outline key approximations, and establish the 
CFTD-SALT correspondence. In Section II, we provide a 
comparative study of CFTD and SALT for a two-mode 
quasi-degenerate laser in two different regimes of param¬ 
eters. Keeping all other parameters the same, we an¬ 
alyze the steady-state dynamics of this laser for small 
7 || = 0.001 ( 7 ||/A = 0.0038) where the SIA is valid, and 
then for yy = 1 (yy/A = 3.8037) where SIA can not be 
guaranteed. Indeed, in the latter regime we illustrate 
that the inversion is non-stationary for a narrow range of 
pump powers, and show how this destabilizes the station¬ 
ary emission and ultimately triggers the synchronization 
of the two modes, to return to a dynamical regime where 
inversion is again stationary. 


I. NON-HERMITIAN SPECTRAL APPROACH 
TO LASER DYNAMICS 

We start with the following form of the Maxwell-Bloch 
equations [2] for the scalar electric field amplitude E{r,t), 
polarization P{r,t) and inversion density D{r,t): 

v2F;+-^ij+ = AioP+, (1) 

P+= -{ina+-f±)P+-i^E+D, (2) 

D = 711 [Do(a:) -D]+ *^[A+(P+)* - (F;+)*P+]. (3) 

Here E = + E~, P = + P~ and we used the 

rotating wave approximation (RWA), valid when the fre¬ 
quencies of the aforementioned fields (~ Da) are much 
larger than their relaxation rates (controlled by yy in 
Class A and B lasers [32] )j typically well-satisfied in the 
optical regime. The laser cavity is characterized by the 
complex-valued refractive index distribution n{r). We 
have in mind a quasi-2D geometry in which case the 
scalar field E{r,t) denotes the z-component of the elec¬ 
tric field for transverse magnetic (TM) polarization, and 
n(r) represents the effective index |35j . In the inversion 
equation, Do{r) represents the possibly spatially inhomo¬ 
geneous pump distribution. We note that the description 
in Eqs. (|2][^ is sufficiently general to describe the salient 
features of various gain media characterized by a single 
dominant optical transition frequency, including quantum 
cascade-based lasers (see Supplementary Information in 
HU]). The remaining parameters are as follows: 7 _l and 
yy are the polarization and inversion decay rates, is 
the center frequency of the gain curve, g is the dipole 
moment of the individual two-level emitters forming the 
gain medium, qiQ is the magnetic permeability, and c is 
the speed of light. 

In the standard spectral approach [3U], the electric 
field and polarization are expanded in a complete set 
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of states e.g. E+{r,t) = Cm(t)(^m(r'), with 

(j)m satisfying V^(j)m{r) = -in?{r) (w^/c^) 4>m{r) with a 
boundary condition at the cavity walls dD that gives rise 
to a Hermitian boundary value problem, and hence a 
complete set of orthogonal states {4>Tn} with real-valued 
frequencies {ujm}- There are two crucial shortcomings 
of this approach. The first is that a phenomenological 
decay rate has to be added by hand to the equations 
m in order for a well-dehned steady-state to exist. As 
an additional consequence, there is no systematic way to 
extend the solution to the exterior of the cavity, where the 
fields are actually measured. A second shortcoming with 
this approach is that spatial hole burning interactions 
can only be captured perturbatively in the electric field 
amplitude, or else through an adiabatic elimination of the 
gain medium degrees of freedom. 

Here, we extend this spectral approach to a consistent 
mathematical framework, by first expanding the electric 
and polarization fields in terms of CF states [18] through 
the following ansatz: 


E+{r,t) ='^enit)(pnir,na)e (4) 

n 


The biorthogonal set of CF states {(pjn{r,^)} is the 
solution to the Laplace eigenvalue problem 'S/'^tprn = 
—n^(r) iprn with outgoing boundary conditions 

~ i{^l/c)(pjn as r- —>• oo. The set of CF states is 
the exact non-Hermitian basis to expand the fields in an 
arbitrary open geometry described by n(r) that is excited 
by an arbitrary spatial distribution of monochromatic 
sources at frequency fl [20] • The solution to this bound¬ 
ary value problem leads to a complex-valued spectrum 
uim and associated eigenmodes (pm that parametrically 
depend on the excitation frequency fl (see Fig. for an 
example of this parametric dependence). The imaginary 
part of LOm provides the crucial mode-dependent losses, 
either through optical leakage out of dD or the material 
absorption described by the imaginary part of n{r). 

A crucial computationally important detail here is that 
the computational domain of the CF problem can be 
reduced to a ’’last scattering surface” dD that can be 
chosen to be the minimal volume that includes all the rel¬ 
evant scattering elements. In practice m, dD is chosen 
to be the minimal circular boundary (in 2D) that includes 
all the spatial inhomogeneities of n(r). Therefore, by 
construction the relevant open boundary conditions are 
exactly satisfied through the use of the CF basis in the 
expansion Eq. Q. In addition, CF states can be analyti¬ 
cally continued straightforwardly outside dD and hence 
the fields, and in particular the electromagnetic flux and 
the measured spectrum, can be calculated exactly in the 
farfield [TO] . 

In the ansatz (0-® the time-dependence of each field 
variable is entirely encapsulated in its respective coeffi¬ 
cients £n and pn- For computational efficiency, we factor 


out the fast oscillation at atomic frequency Da- The 
spatial dependence is entirely captured by the CF states, 
which are calculated, in a departure from previous applica¬ 
tions of the CF basis, only at Dq. This is a very good and 
well-controlled approximation, for the CF states and fre¬ 
quencies {prn{'t', D),a;m(D)} typically change slowly when 
the excitation frequency D is varied, see for an example 
Fig-iD This is in fact one of the crucial factors in the 
computational efficiency of SALT m- We note that it 
is possible to choose an unusual geometry where for a 
certain narrow regime of parameters (D) this assumption 
may fail, but generally this should be taken as a hint 
that some extraordinary spatial physics is present in the 
system that may give rise e.g. to an exceptional point 

m- 
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FIG. 1. (a) CF Eigenvalues and their variation as a function 
of the external frequency are shown for a ID cavity with 
n = 3, Da = 20, 7 _l = 1, and nQaLjc = 60 (see description of 
parameters in text). Closely spaced blue markers indicate 
the variation in eigenvalues within the range D = Da ± 7x 
and the red ’x’ marker indicates the eigenvalue calculated 
at D = Da. The shift in the real part of these eigenvalues 
is negligibly small, (b) The variation in the CF eigenstate 
associated with the smallest eigenvalue shown in (a). Blue lines 
show the eigenstate calculated in the same range as (a) while 
the red line shows the eigenstate calculated at D = Da. The 
inset zooms in on the peak enclosed in the dotted rectangle, 
closely showing the very small variation in the eigenstates as 
a function of D. 


With the above ansatz of Eqs. (|4][^ inserted into Eqs. ([^ 
1^, we can derive the following equations of motion for 
the time-dependent dimensionless coefficients e(t), p(t), 
Djyin (t') . 




Pm 


b 


mn 




n„ 


mn Pn 5 


- 7-LPm - *7-L X! ^ 

n 

7||(-Do ,mn Emn) 

2 


+ * „ ( ^ '^ ^mrsnSrPs ^mrsn^rPs)' 


( 6 ) 

(7) 

( 8 ) 


Here we introduced the inversion matrix Dmn{t) = 
/cavity 'n{r)p^{r)D{r,t)pn{r), a set of space- 

independent coefficients describing the mode-projected 
inversion distribution. While the inversion D{r,t) itself is 
real-valued, the coefficients are in general complex¬ 
valued. All the variables are rendered dimensionless 
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tlirOUgh Sm — ^mj^ct Pm — PmjPct D^yi^i — DYnnj 

using the following scale factors that contain all the units: 


Er_ = 


25 


Pr = 


Ec 


Dr = 




(9) 


Furthermore, time and decay rates are scaled by the effec¬ 
tive cavity round-trip time tuT = nLjc (n can be taken 
to be the spatially averaged effective index, and L = 
for a cavity with volume V). The key step in obtaining 
Eqs. (|6][^ is the elimination of the spatial dependence of 
each field vector by utilizing the biorthogonality of the 
CF basis vectors. We will drop the tildes henceforth. 


/ cavity 


dr n{r)tpm{r)ipn{r) = 


( 10 ) 


In contrast to a Hermitian orthogonality relation, this 
inner product does not contain a complex conjugation. 
This is a consequence of the dual modes (left eigenvectors) 
{(pm} satisfying the relationship (pm = pffn. [ISlIini- 
This step produces the following unitless complex¬ 
valued parameters appearing in the above equations. 


'^mrsn 


= L 


/cavity 


dr n{r)^pm(r)ifr(r)(p*ffr)pri{r), 


( 11 ) 


A' = T. 

■^mrsn ^ 


dr n{r)(pm{r)(p*(r)(pffr)p>n{r), 


' cavity 


^mn — 

f 


/cavity 

^0,mn — 

/ 


/cavity 

Here, Amrsn 

and A‘, 


dr ipm{r)(pn{r), 


dr n{r)(pm{r)Do{r)(pn{r). 


( 12 ) 

(13) 

(14) 


can be seen as a generalization 
of the inverse mode volume in the Hermitian version 
of the single-mode laser problem. Interestingly, Bmn 
is not diagonal unless the index is uniform across the 
cavity. The effective mode-projected pump parameter is 
given by Flo.mn and is the most critical parameter here. 
These overlap integrals Eqs. (TTpi) are calculated prior to 


numerically solving the time-dependent system of coupled 


equations in Eqs. (6 8) and they encapsulate the impact 


of the resonator modal structure on modal interactions. 

An important aspect of the above spectral formulation 
of semiclassical laser equations is that it takes into ac¬ 
count modal interactions through spatial hole burning 
exactly. Majority of the past spectral methods (with the 
exception of SALT), account for interactions only pertur- 
batively and generally to third order in the electric field 
amplitude. This approximation, as pointed out first in 
[55] and later in quantitative detail discussed in |^, is 
only valid near the lowest laser threshold, and generally 
severely underestimates the number of lasing modes at 
higher pump powers. 

As long as the parametric variations of the CF basis 
is small within a window 'y± of fla, the Eqs. (6][8) are 


exact up to the slowly varying envelope approximation 
used in Eq. to remove second order time derivatives 
in in and The impact of the latter in SALT has been 
quantified perviously m and was shown to introduce 
small inaccuracies in the calculation of steady-state lasing 
characteristics, but was not found to lead to any qualita¬ 
tive differences. This approximation is not critical to the 
success of the method as discussed below, and can easily 
be undone at the expense of introducing additional fields. 

In the next section, our goal is twofold. We would 
first like to benchmark SALT against the CF-projected 
time-dependent laser equations Eqs. ([6][^ (CFTD) in the 
regime of parameters where SALT is known to be accurate. 
Next, we investigate a regime accessed by the change of 
a single parameter, yy, leaving all other parameters the 
same, where the SI A is suspect. Here we encounter a 
narrow regime of pump powers where the system is critical 
and unstable towards a synchronized oscillation regime. 
In this regime that, for the special cavity configuration 
of Fig. occurs at extremely high pump power (about 
25 times the lowest threshold), SALT fails to capture the 
underlying dynamics qualitatively. Interestingly, below 
and above this narrow regime of pump powers, the SIA 
is valid and SALT is accurate. 

A second aim of the following discussion is to present an 
accurate picture of the synchronization transition, known 
as cooperative frequency locking [55]. Our theoretical 
result is able to accurately capture the interesting dynam¬ 
ical regime around the critical pump power for locking, 
experimentally observed for the first time in 1988 [39] . 


II. BENCHMARKING SALT AGAINST CFTD: 
THE TWO-MODE QUASI-DEGENERATE LASER 


In this section, laser dynamics is investigated for a quasi¬ 
degenerate ID cavity. It consists of a dielectric slab with 
refractive index n = 3.3 which sandwiches symmetrically 
a layer of index n = 1.5 of thickness SL (see Fig. [^. 
The gain curve is centered at fla = 20.5, 7 _l = 8, and 
nVlaLlc « 135. This particular choice of 7 _l ensures a 
flat gain experienced by both of the cavity resonances 
included in the calculations below, significantly reducing 
the effect of gain-pulling in the time-dynamical scenario 
where lasing mode frequencies shift strongly with the 
pump power. 

Choosing SL/L = 1/40, the resonances of the cavity 
come in quasi-doublets which are separated from each 
other by a relatively large spectral range (See Fig. [^. 
These conditions are ideal to consider two regimes, one in 
which the SIA is valid (Regime A) and another where it 
cannot be guaranteed (Regime B), by changing the value 
of a single parameter, yy, and leaving all other parameters 
identical. 

In regime A (yy ^ A ^ 7 _l), the assumptions under¬ 
lying SIA are valid and SALT and CFTD results should 
agree quantitatively [57|. We will first set up the cor¬ 
respondence between SALT and CFTD variables in the 
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FIG. 2. (a) The refractive index distribution of the cavity 
described in the text is shown as a function of space (r), 
scaled by cavity length L. The CF eigenvectors color-matched 
with the CF eigenvalues marked in (b) are plotted inside and 
outside the cavity; dashed gray lines on each end mark the 
cavity faces. Note that the CF mode with an anti-node in the 
low refractive index region shows the expected large magnitude 
due to reduced index, (b) The 10 eigenvalues closest to gain 
center are shown for the cavity shown in (a); filled red/blue 
circles mark the two eigenvalues used in the calculations in this 
section: uji = 20.3764 — 0.0970i and uJ 2 = 20.6393 — 0.0957i. 
The spacing between them is A = Re [tU 2 — tui] = 0.2629 and 
their spacing from the adjacent eigenvalue pairs is A 2 « 0.7. 


steady-state and then demonstrate excellent agreement be¬ 
tween the two methods using the 2-mode quasi-degenerate 
laser as an example. 

In regime B (A ^ yy 7 _l) however, accessed here 
by changing yy, the stationary inversion approximation 
can not be guaranteed. Indeed, while at low powers the 
laser oscillates in two frequencies (“two-mode lasing”), 
above a critical pump power Dq = corresponding 

to the threshold for synchronization^ these two frequencies 
lock and a single frequency remains. Just prior to syn¬ 
chronization, a close up at the power spectrum of various 
dynamical variables of CFTD in Eqs. (6]|8) reveals that 
close to the synchronization threshold the SIA breaks 
down. The SIA remains valid generally however, breaking 
down in a very interesting way but only within a narrow 
range of pump powers. 


SALT-CFTD correspondence 

The strength of a steady-state approach like SALT [18] 
is that it directly delivers the frequencies as well as the 
intra- and extra-cavity held amplitudes as functions of the 
pump power Dq. As such, it is not immediately clear how 
the SALT variables are related to the CFTD variables 
£m, Pm and Dmn- In this subsection, we will set up 
this correspondence when this correspondence exists, and 
then compare SALT and CFTD results in the following 
subsection. 

SALT is obtained by making a more specihc ansatz 
for the long-time solution of MB equations than that for 


CFTD: 


£;+(r,t) = ^T('^)(r)e-*^'''’‘, 

l-L 

(15) 

P+(r,t) = ^p(^)(r)e-*^'''’*. 

(16) 




The crucial point here is the assumption of a specihc form 
for the exact time-dependence once steady-state is reached 
(compare Eq. (15) to Eq. (|^). The helds are assumed 
to be expandable in a discrete Fourier representation 
with a hnite number of laser frequencies which are 
unknown and to be determined. 'I'(^^(r) are the spatial 
held amplitudes corresponding to the exact (non-linear) 
lasing modes, also to be determined through the SALT 
equations: 

+ eg{f))nl] ^ f,{r) = 0, re cavity (17) 
7-L Doir) 


= 




(18) 


Here F^ = Ty/bl + — Ha)^]- Note that the polariza¬ 

tion spatial amplitudes (r) can be directly related to 
fl'('^)(r) and do not show up in the hnal set of equations 
to be solved. Also note that in contrast to CFTD, the 
index /i specihcally identihes lasing modes oscillating at 
distinct frequencies (as opposed to spatial modes). 
The time-independent SALT equations Eq. (17) are then 
solved by projecting each lasing mode (r) into a set 
of CF states for the associated frequency of oscillation 


(19) 

n 

The SALT-CFTD correspondence is unveiled by assuming 
(/3ji(r, « (/?„(r,Ha), which as discussed before, is 

generally a good approximation. In that case, 

£^{t) = e*^“‘ (20) 


where the coefficients on the left and right hand side are 
the CFTD and SALT variables, respectively. 

Additional insight is obtained by asking what assump¬ 
tions SALT makes about the solution of CFTD equations 
Eqs. (6]|8), that are more general. SALT corresponds to 
specific long-time solutions of the CFTD equations for 
which e„(t) = Eu Pn{t) = and 


Dmn = 0. The last assumption is the mode-projected 
version of the SIA and one of the consequences is that yy 
drops out of the equations. That doesn’t however mean 
that SALT solutions do not depend on yy, but rather that 
the entire yy-dependence of SALT solutions is contained 
in the particular scaling of the electric field Eq. Of 
course, being exactly equivalent to MBE equations up to 
the aforementioned approximations, the CFTD equations 
permit far more general solutions, one of which we will 
encounter further below. 
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For an A^-mode CFTD calculation where N is the num¬ 
ber of modes, a CF basis of equivalent dimension must 
be constructed, and -|- 2N equations must be solved. 
For a 2-mode calculation, this amounts to 2 equations 
each for the electric and polarization fields, and a total 
of 4 equations for the diagonal and off-diagonal elements 
of inversion. Below we will discuss the two-mode regime 
for the quasi-degenerate ID laser we introduced before 
(Fig.[2|. 

Regime A: Stationary inversion 

We use the parameters quoted at the beginning of 
Section II and take qy = 0.001. Right at the onset of the 
second mode this gives qy/A = 0.0038. A only slightly 
changes in the calculated interval of pump powers (see 
Fig.j^b)) and the assumptions underlying the SIA remain 
rigorously valid throughout. In the figures below, we use 
a normalized pump (3 = Dq /Dth where Dth is the lasing 
threshold. 

CFTD calculations show that both |£i(t)p and |e 2 (t)P 
reach steady state after initial transients die out. Some 
sample time-series are shown in Fig. |^a,b) for two dif¬ 
ferent pump powers [3 = 1.16 and (3 = 2. These findings 
indicate that there is a single frequency in both £i (t) and 
£ 2 ( 0 , as conhrmed in the retentive power spectra shown 
in Fig. [sj^c). Shown in Fig. [^d), \Djnn{^)\‘^ is the power 
spectral density (PSD) of Dmn{t) (itself not shown). Here 
we plot |ZIii(a;)p and |Di 2 (a;)p only. Further detail is 
shown in the inset which zooms out and shows in loga¬ 
rithmic scale that the non-stationary components in D 12 
(red peaks) are suppressed by more than three orders of 
magnitude with respect to the static component of Du 
(black peak). The smallness of these side-peaks indicates 
that the SIA is an excellent approximation and the SALT- 
CFTD correspondence should be possible, which is what 
we do next. 

The SALT calculation containing two lasing modes 
expanded into a basis of two CF eigenvectors will con¬ 
tain four coefficients (a$^^\ a^\ o!'i \ 0 -^ 2 '^) and two lasing 
frequencies As discussed in the previous 

section, in the steady-state, the information contained 
in these SALT variables can be retrieved from the two 
time-dependent CFTD variables (£ 1 , £ 2 ). To do so, we 
simply expand and rearrange the SALT ansatz for two 
lasing modes, 

E+{r,t) = (a^e-*^'"'* +al')e-*^'''‘)(^i(r,L!,) (21) 

-I- -I- 4^^e“*^'^’‘)(p2(^, ^a) (22) 

= Ha) -h £2it)(p2ir, Ha) (23) 

The CFTD results imply that £„(t) « in 

other words « 0 for n ^ m. In SALT language, 
this means that the single-pole approximation is valid 
throughout the calculated regime - the CF eigenvectors 
calculated for a cold cavity very closely represent the two 


lasing modes g single CF component is 

sufficient to represent each mode. For CFTD-SALT com¬ 
parison, in Fig. Qa) we plot the “intensity”, J2n 

from SALT, and compare it to Im = T 2 -T 1 Iti km(^)Pi 
calculated for a sufficiently long sampling time after the 
steady state is reached in CFTD. 



Frequency (a) Frequency (o) - 


FIG. 3. (a) The time-domain behavior of |£i(t)p (blue) at 
(3 = 1.16 in the single-mode lasing regime, (b) The time- 
domain behavior of |£i(t)p (blue) and |£ 2 (t)|^ (red) in the 
multi-mode regime. Their approximate steady-state values 
are labeled, (c) The PSDs of £\{t) (blue) and £ 2 )^) (red) at 
the same (3 values as above; the labeled values compare well 
to those marked in (a) and (b) and in Fig. ga). (d) The 
spectral content of the diagonal inversion element Du (black 
peak); inset log-plot also plots the off-diagonal element D 12 
(smaller red side-peaks) and shows that it’s nearly 3 orders of 
magnitude smaller than Du. 

The threshold of the first mode as calculated by SALT 
and our time-dynamical method is almost the same, and 
the emission intensities also coincide up to the point where 
a second mode begins to lase in the SALT calculation [see 
Fig. B a)]. Shortly thereafter, the second mode begins 
to lase in the time-dynamical calculation as well and 
both modes progress with comparable slope-efficiencies 
up to high pump powers. The steady-state frequencies 
[Fig. I^b)] confirm the expected steady-state behavior. 
The small offset between the CFTD and SALT solutions 
can be attributed to the use of SVE in CFTD, whereas 
SALT does not make this approximation (See Ref. [37] 
for a discussion of this point). 

Regime B: Non-stationary inversion 

Keeping all other parameters, we now choose yy = 1. 
At the onset of the second mode, this gives yy/A = 5.61. 
We will find that A will change dramatically in this case, 
essentially going to zero as the pump power is increased. 
This phenomenon is known as ’cooperative frequency 
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FIG. 4. (a) The steady-state emission intensity of the two 
lasing modes as calculated by SALT (red) and CFTD (blue). 
The values shown here are a time-average (details in text) of 
the steady-state time domain behavior at each pump step; 
the time-domain behavior is shown in Fig. [^a-b) for selected 
values of /3 = 1.16 and (3 = 2 (marked in this plot), (b) The 
steady-state center frequency for both lasing modes for SALT 
(red) and CFTD (blue). 


locking’ [29], and has been experimentally studied in 

Ref. [55] . 

In Fig.[^ CFTD reveals that while the first mode starts 
lasing at the same threshold as before, with a nominally 
identical lasing frequency = 20.3762, the second mode 
lases at a threshold nearly 10 times larger with frequency 
O 2 = 20.5426. However, immediately after the turn-on 
of the second mode, |£:„(t)p ceases to reach a stationary 
value, implying the existence of multiple frequencies in the 
respective spectra e„(a;). In lieu of intensities, we plot here 
the time-averaged quantities Im, and indicate the size of 
oscillations, 5, around the mean by shaded regions. As the 
pump approaches the synchronization threshold 
the oscillations in the intensities grow (for the second 
mode, the oscillation magnitude remains always of the 
order of the mean, implying a clear limit cycle solution). 
For Dq > the oscillations in the intensities abruptly 

disappear, and all field amplitudes oscillate at a single, 
synchronized frequency. The synchronization threshold is 
clearly defined and corresponds to (3 = 22.98. 

A close look at the dynamical behavior of the system 
near the synchronization threshold is provided in Fig. [^ 
We follow here the PSDs of (a) the electric field and (b) 
the Dll and D 12 components of inversion as a function of 
the pump. Note that these spectra are shown for only a 
small range of pump powers around the synchronization 
threshold at /3 = 22.98. The largest peak in Fig, j^a) 
belongs to the dominant mode, shown in red in Fig. ^a), 
and the sub-dominant peak belongs to the mode shown 
in blue in Fig. [ga). As pump power is increased, A 
decreases, and additional peaks enter the monitored fre¬ 
quency window, separated by integer multiples of A (with 
respect to the original laser frequencies fl^). These are the 
frequency-doubling harmonics mentioned in |39j . Note 
that the highly non-linear sawtooth-like oscillations in the 
intensities seen in Fig. [^b) are closely linked with this 
proliferation of frequencies in the power spectrum. As the 
pump power is increased further, all these peaks approach 
each other in a dramatic manner and at /? = 22.64, recol- 


FIG. 5. (a) The steady-state emission intensity of the two 
lasing modes. The first mode (red) reaches threshold at the 
same pump power as Fig. while the second mode (blue) 
reaches threshold at more than ten times the threshold of the 
first mode. At P = 21.91 the two modes begin to converge, 
leading to a sharp rise (decline) in the blue (red) mode ac¬ 
companied by oscillations steadily increasing in amplitude 5 
(shown by the light red and blue shaded triangular regions). 
After synchronizing, the modal intensities of the two modes 
are comparable and they increase linearly, (b) The top figure 
shows the large oscillations in time in each mode, immedi¬ 
ately before synchronization; the straight lines identify the 
time-average of the oscillating signal and correspond to the 
values marked in (a). The bottom figure shows the absolute 
absence of these oscillations and a return to a steady-state 
after synchronization. 

feet into the single peak shown (in purple) at (3 = 22.98 
and beyond. This peak is seen to be shifted from the 
point of convergence and from both primary frequency 
components, and it is pulled towards the gain center. A 
slightly different perspective is offered by the evolution of 
the power spectrum of the inversion [Fig.|^b)] which also 
shows that the off-diagonal frequency components (blue) 
converge into the DC component (red) as they must if 
there is to remain only one mode. The new mode that 
emerges beyond synchronization is comprised of nearly 
equal contributions from both CF states, which can be 
seen directly from the coming together of |eip and |e 2 p 
in Fig. I^a). The new mode has a non-trivial spatial pat¬ 
tern, which is embodied in a non-linearly generated phase 
between the two CF states composing the new laser mode. 
More detail on this point is provided in the Appendix. 

It is interesting to see how all this looks from the 
perspective of SALT. We provide a comparative study in 
Fig. HI The SIA appears to be valid everywhere outside 
the comparatively narrow range of pump powers 15 ^ 
(3 < 22.64 and the SALT-CFTD correspondence should 
in principle be possible. Comparing the intensities in 
Fig. 0 (a), however, we see what is a strinkingly large 
discrepancy between SALT and CFTD for (3 « 15. While 
SALT finds two modes that turn on relatively close to 
each other {D^ « l.SD^^), CFTD shows that the second 
mode does not turn on before (3 = 12.86. These seemingly 
disparate results should however be taken with a grain 
of salt. We first point out that the two thresholds found 
by SALT are identical to those found for yy = 0.001 
shown in Fig. This is of course expected because SALT 
equations do not depend on yy when expressed in scaled 
variables (Eq. which is what is plotted in the vertical 



















FIG. 6. (a) All red (blue) peaks belong to the spectrum 

of El (£ 2 ). All peaks are seen to draw closer and converge 
at /3 = 22.64 (marked by a star) and a single purple peak 
representing the nearly identical spectra of ei and £2 can be 
seen at /3 = 22.98 = Dly„i. (marked by a pentagon) and beyond, 
(b) The red (blue) peaks belong to the Dn (O 12 ) component 
of inversion. Compared to the inset in Fig. I^d), the peaks of 
Di 2 are comparable in magnitude to the DC peak from Du, 
representing the significant effect of time-dynamical behavior 
in this calculation. Similar to (a), all peaks converge upon 
synchronization and only a DC component (purple) remains 
from Dll and D 12 . 


axis. The large discrepancy (despite SIA appearing to 
be valid) is simply because SALT predicts the turn-on 
of the second mode incorrectly, by a large margin. The 
consequence is that the change in slope of the intensity 
of the first mode that happens when the second mode 
turns on is incorrectly predicted by SALT as well. The 
seemingly large discrepancy between intensities by the 
time the second mode turns on in CFTD at /3 = 12.86 
is thus simply due to the incorrect slope. We note that 
the pump range we are comparing is extremely large (and 
could well be inaccessibly large for certain gain media) 
- the synchronization threshold found is about 22 times 
larger than the (lowest) laser threshold. 

The culprit for the incorrect prediction by SALT of 
the threshold of the second mode is interestingly still due 
to the breakdown of SIA, but in a non-trivial manner. 
While the oscillating corrections to the inversion are still 
small for /3 < 15, they generate a polarization component 
oscillating at ili + A i.e. 122, that is proportional to the 
intensity of the first mode ~ |'I'*^^)(r)p that does become 
large as pump power is increased. It can be shown that the 
threshold condition of the second mode is changed by a 
term proportional to ( 7 ||/A)|'k(^^(r)p. An appropriately 
modified set of two-mode SALT equations can be found 
and its implementation would correctly reproduce 
the behavior seen in CFTD for /I < 15. 


However, SALT will have nothing to say and will fail 
qualitatively in capturing the physics in the range of pump 
powers plotted in Fig. very near the synchronization 
threshold. This is directly linked with the appearance of 
oscillating terms in the inversion (see Fig. |^b)) that are 
comparable in magnitude to the static terms. Note that 
very interestingly the oscillations only appear in the off- 
diagonal elements, while diagonal elements mostly remain 
stationary. An analytic understanding of these features 
will be investigated in future work. 




5 10 15 20 26 30 34 38 


Pump (P) 


FIG. 7. SALT and CFTD calculations are compared for 
the same 2-mode calculation as above both before and after 
synchronization (yy = 1). (a) A comparison of the intensities 
as a function of pump. CFTD calculation for intensities of 
the two lasing modes is shown in red/blue circles while the 
SALT result is shown in red/blue crosses. Purple triangles 
represent the sum of the intensities of the two modes in the 
time-dependent calculation. The shaded green (white) region 
marks the pump range before (after) synchronization. Black 
squares show the SALT calculation after synchronization, (b) 
A comparison of the frequencies as a function of pump. Again, 
red/blue circles mark the time-dependent calculation and 
red/blue crosses mark the SALT calculation. Purple triangles 
(black squares) show the time-dependent (SALT) synchronized 
solution. The two smaller plots zoom in on regions of interest. 

Post-synchronization, as seen in Fig. for j3 > 
a properly conditioned SALT (discussed m the Appendix) 
very precisely predicts the synchronized mode, both its 
oscillation frequency and the spatial composition. This 
again, is not surprising because now the SIA is valid to an 
excellent approximation in a single-frequency regime of 
lasing. We note however that SALT is unable to capture 
the synchronization threshold accurately. 
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III. CONCLUSION 

Here we’ve presented a computationally and concep¬ 
tually efficient approach to isolating and studying time- 
dependent effects in lasers. Using a spectral approach, 
we fully treat the open nature of lasers and integrate 
out the spatial variables, obtaining dynamical equations 
for the time-dependent coefficients describing the electric 
and polarization fields and the inversion. This delivers a 
highly scalable multi-mode framework for analyzing in¬ 
trinsically non-stationary phenomena in open resonators 
of arbitrary spatial complexity, gain medium distribution, 
and pump prohle. The simplest of such effects, mode 
synchronization, is studied here in a simple ID cavity 
featuring pairs of closely spaced quasi-degenerate modes 
(small A). With small enough yy, we obtain a stationary 
behavior once the transients die out, as postulated at the 
outset. At larger yy, non-stationary behavior is demon¬ 
strated in a narrow range of pump powers. We find that 
the stationary inversion approximation is largely valid 
in the parameter regimes investigated here, failing only 
in a narrow range of pump powers, for very large yy/A, 
and for a special choice of the resonator structure. We 
expect that CFTD will find application in particular in 
modeling time-dependent phenomena in quasi-2D and 3D 
laser structures, because of its efficient spectral decom¬ 
position method that takes into account the openness of 
the underlying resonator structure in essentially an exact 
manner. 
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V. APPENDIX 

In this Appendix, our goal is to provide more detail on 
the SALT-CFTD correspondence in Regime B. The SIA is 
valid to a good approximation for /3 < 15 and (3 > 22.98, 
and a SALT-CFTD correspondence in these power ranges 
is therefore possible. 

As discussed in Regime B and Fig. above, in the pre¬ 
synchronization regime the apparent sizable discrepancy 
between SALT and CFTD solution is understood, and 
can be accounted for by a modified version of SALT m- 
We focus here on the SALT-CFTD correspondence in the 
synchronized regime, where it is important to properly 
condition SALT. 

The standard SALT algorithm uses an ‘adiabatic’ sweep 
of pump power (not in the dynamical but computational 
sense). In other words, the solution in the previous step of 
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FIG. 8. (a) Frequencies as a function of the pump are shown 
for the seed (black triangles), and the three SALT solntions: 
the single-mode solution near the cavity frequency (blue cir¬ 
cles), the synchronized solution (green diamonds), and the 
unstable solution (red squares), (b) Intensities are shown 
as a function of pump for the same color-scheme as (a), (c) 
and (e) The distribution of seeds in phase space is shown for 
£i and £ 2 ; inset in (c) zooms in a region near the real axis 
where the synchronized and unstable solutions are found, (d) 
SALT solutions for £ 1 , which lie on the real line due to the 
SALT gauge condition, (f) SALT solutions for £2 and their 
distribntion in phase space. 


pump power Dq is fed as a seed for the non-linear solver 
for the next pump power. This practical procedure speeds 
up the computation in a dramatic way. If this were done 
blindly, then the two SALT solutions found in Fig. for 
(3 < would simply extend without any apparent 

discontinuity to larger values of pump power, missing the 
synchronized solution. In fact, SALT has multiple single¬ 
frequency fixed points for [3 > two of these are 

composed of a single CF mode (i.e. single pole) and the 
other two are composed of a particular balanced superposi¬ 
tion of two CF modes (multi-pole solution). Interestingly, 
one of the latter is the synchronized solution found in the 
long-time limit of CFTD (shown using green markers in 
Fig.§. This indicates that SALT does capture the fact 
that there is a single stable oscillation frequency and that 
the spatial structure of this laser mode is such that it is 
a particular superposition of the two spatial modes that 
were oscillating independently at lower powers. Fig.[^d,f) 
indicate that three of these fixed points are stable, and 
one is unstable, as revealed with different initializations 
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of the SALT non-linear solver (Fig. |^c,e)). The unstable 
solution is a synchronized solution that is orthogonal to 
the stable synchronized solution. The stability of the two 
’’single-pole” solutions (only one shown in the frequency 
window plotted) from the point of view of SALT is a per¬ 


ceived one, and is due to the neglect of the non-stationary 
terms in the inversion that in turn changes the stability 
structure of the solutions. We conclude that care must be 
exercised when conditioning SALT solutions for regimes 
outside its stated validity, even when the SIA appears to 
be a good approximation. 
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